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ABSTRACT 

We have investigated the effect of an assembly history on the evolution of galactic dark matter 
(DM) halos of 10 1 h~ x MQ using Constrained Realizations of random Gaussian fields. Five different 
realizations of a DM halo with distinct merging histories were constructed and have been evolved using 
collisionlcss high-resolution Y-body simulations. Our main results arc: A halo evolves via a sequence 
of quiescent phases of a slow mass accretion intermitted by violent episodes of major mergers. In the 
quiescent phases, the density is well fitted by an NFW profile, the inner scale radius R s and the mass 
enclosed within it remain constant, and the virial radius (i? v ir) grows linearly with the expansion 
parameter a. Within each quiescent phase the concentration parameter (c) scales as a, and the mass 
accretion history (M vu .) is well described by the Tasitsiomi et al. fitting formula. In the violent phases 
the halos are not in a virial dynamical equilibrium and both R s and i? v ; r grow discontinuously. The 
violent episodes drive the halos from one NFW dynamical equilibrium to another. The final structure 
of a halo, including c, depends on the degree of violence of the major mergers and the number of 
violent events. Next, we find a distinct difference between the behavior of various NFW parameters 
taken as averages over an ensemble of halos and those of individual halos. Moreover, the simple scaling 
relations c— M v ; r do not apply to the entire evolution of individual halos, and so is the common notion 
that late forming halos are less concentrated than early forming ones. The entire evolution of the halo 
cannot be fitted by single analytical expressions. 

Subject headings: cosmology: dark matter — galaxies: evolution — galaxies: formation — galaxies: 
halos — galaxies: interactions — galaxies: kinematics and dynamics 



1. INTRODUCTION 

The hierarchical buildup of cold dark matter (CDM) 
halos is a strongly nonlinear process. It is associated 
with the buildup of a unique density profile — one of the 
most fundamental charac teristics of the DM ha los. Based 
on Y-body simulations, iNavarro et alJ l)1997l hereafter 
NFW) found that the CDM halos can be universally fit- 
ted by a two-parameter functional form: 

o( r ) = — (1) 

P[r> (r/Rjil + r/R,)*' [) 

where R s is a characteristic "inner" radius at which the 
logarithmic density slope is —2 and p s is the density at 
R s - The cosmological evolution of the NFW parameters, 
therefore, is an issue of a broad interest and is addressed 
in this work. 

A useful alternative parameter to describe the shape 
of the profile is the so-called halo concentration pa- 
rameter defined as c = i? v i r /i? s , where i? v ir is the 
halo virial radius. Although such a profile has been 
confirmed by numerous Y-body simulations, the ex- 
act value of the inner slope parameter is uncertain. 
Several authors hav e found that halos have d ensity 
cusps steeper (e.g.. iFukushige fc Makinol 119971 1200.4 
iMoore et al.lll999HGhigna et al.ll2000H or shallower (e.g., 
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iSubramanian et al J 12005 iTavlor fc Navarro! I2001|) than 
-1, the NFW value. The NFW profile and its modifica- 
tions are specifi c cases of a three -parameter profile fam- 
ily proposed_by^femfluis3 ||1990), an d further developed 
bvlZhaol (I1996IK On the o ther hand. llJing fc Sutoll2000l 
l2002j) . lKlvDin et all lj200l and lRicotti l |2003ft found that 
the CDM halos do not maintain the universal density 
profiles, and the inner slope chan ges from galaxy-size ha - 
los to cluster-size halos (see also iTasitsiomi et alJl2003l . 

These numerical results appear in conflict with the 
observational evidence — rotation curves of low sur- 
face brightness galaxies yield d ensity profiles with nearly 
constant density cores (e.g., iFlores fc Primackl 119941: 
iSalucci fc Burkertl I2000J) . St udies of brighter galax- 
ies imply similar problems llSalucci fc BurkertJ 120001 
Ide Blok fc Bosmal2002HGentile et all2004j) . While some 
studies of the weak l ensing seem to suppo rt the NFW 
density profile (e.g., iHoekstra et al.l l2004|) . in a spe- 
cific case of the cluster A1689, the required concentra- 
tion parameter is dramatically larger than the typical 
value obtained in the si mulation of a cluster-size object 
l)Broadhurst et al]l20l)5|) . 

In principle, one can list three different options to 
resolve the above difficulties. First, the inner few-kpc 
rotation curves of disk galaxies can be contaminated 
with the (unresolved) won -circular motions triggered 
e.g., by stell ar bar s (e.g.. iBlais-Ouellette et al J 120011 
i Bolatto et all l200l ISimon et all 120031 iWeldrake et alJ 
120031 iSimon et alJl2005|) . Second, the NFW approach 
is to neglect the triaxial nature of DM halos by "spheri- 
calizing" while analyzing them. Furthermore, the issue of 
the "cusps" might be related to the effective resolution 
in Y-body simulations or other numerical effects (e.g., 



2 



Romano-Diaz et al. 



IMoore et alJll99fift . 

Finally, the difference between the collisionless simu- 
lations and observations can underscore a physical ef- 
fect, like the absence of dissipation in the pure DM mod- 
els and the presence of dissipative baryons in the real 
galaxies. This apparent discrepancy can be reconciled 
within the context of the CDM model by considering the 
effect of clumpy baryons erasing the cusps on relevant 
spatial scales, from g alactic halos to clusters of galaxies 
ijEl-Zant et all I2001L 1200^) . Alternatively, it has been 
suggested th at the stellar bars facilit ate the destruction 
of DM cusp s llWeinberg fc Kat3l20f)2t) . but this has been 
disputed bv lDehnenl l|2005l) . 

There has been no "natural" explanation for the ori- 
gin and "universality" of the NFW profile. While an- 
alytical models necessarily invoke spherical symmetry, 
numeric al simulations emphas ize the background asym- 
metry. iGunn fc Gottl l)1972|) proposed an analytical 
model which invokes the collapse of uniform spherical 
perturbations of collisionless DM in an expanding uni- 
verse. This model explains some of the global prop- 
erties of virialized halos (e.g., mean density and size), 
but does not account for the density profile. A more 
rigorous and exact analytical solution relevant to the 
problem is that of a single scale free spherical density 
perturbation in a Fri edmann universe, the so-called sec- 
ondary infall model (|Gurm1 119771 iFillmore fc Goldreichl 
IT98llBertschingerll985li' and its applica t ion to cosmolog- 
ical models (| Hoffman fc Shah amfll985l: iRvden fc Gunnl 
IT987tlZa,ronbi fc Hoffmanlll 993lLokas fc Hoffmai^OnCfl. 
An orthogonal approach has been suggested to explain 
the emergence of the NFW density profile as the out- 
come of mer gers between substructures and progeni- 
tors of halos llSver fc WhitdllQM INusser fc Shethlll999l 
iSubramanian et all 1200(H) . lEl-Zantl l)2005(l has shown 
that the NFW density profiles remain invariant under in- 
teractions with the subhalos — a necessary step toward 
the universality of this m ass distribution. 

Several authors (e.g., iNusserl [20011 lAscasibar et alJ 
120041 IHiotelisll2002j) extended the secondary infall model 
by allowing for non-radial orbits. They have shown 
that by properly tuning the anisotropy of the orbits the 
model yields a cuspy density profile, quite simi lar to the 
NFW one - a task performed by major mergers l)Lu et alJ 

At least in spherically-symmetric models, it was ar- 
gued that the halo concentration c increases with lower 
mass because they form early, when the uni yerse den- 
sity i s higher. This leads to higher p s (e.g., lEke et all 
1200 1|) . On the other hand, within the hierarchical struc- 
ture formation, halos of a particular mass might form at 
different times, so their mass is not a unique function 
of their formation time, but can depend on their envi- 
ronment. Hence, the central densities do not necessarily 
reflect the properties of the universe at spec ific times 
ijAvila-Reese et all 120051 iWechsler et al.ll2005|) . but the 
si gnificance of this effec t is not clear yet. 

IWechsler et alJ l)2002l hereafter, W02) and others have 
provided simple analytical fits to the evolution of various 
quantities which characterize the halos, e.g., c, halo for- 
mation time, its mass, R s , etc. Yet, the scatter around 
these fits is considerable and its origin is unclear. This 
has led us to embark on a series of numerical experiments 
carefully designed to investigate the evolution of the basic 



halo parameters. By usin g constrained realizations (CRs; 
ijHoffman fc R.ibaklll991ll ) of the initial density field we 
can 'design' the merging history of a single galactic halo. 

Most of the studies of the cosmological evolution 
of DM halos have focused analyzing large ensembles 
of simulated halos and studyi ng their ensemble av- 
erage d properties (e . g., W02 iLemson fc Kauffmannl 
19991 iBullock et al] EiOOlalbl: iPeirani et all 120041 
Avila-Reese et al.l 120051 fReed et al.l 120054 IShaw eHf] 
2005:). Closer inspection of individual halos shows that 
their evolution is not smooth and monotonic but goes 
through a number of discontinuities. The question that 
arises is whether this erratic behavior has a numerical 
origin (e.g., numerical resolution and coarse time sam- 
pling) or is a consequence of the underlying physics. 
In the latter case, this implies that the smooth fitting 
formulae do not provide a good approximation to the 
dynamics and evolution of individual halos. This issue 
is addressed here. 

In a previous work ijRomano-Diaz et al] l|200fiD . here- 
after Paper I) , we have used the CRs of a single galactic 
halo to show that its evolution can be described by a se- 
ries of step functions. By changing the merging history 
in the initial conditions of the same single galactic halo, 
we have found that the inner structure (i.e., within R s ) 
remains unchanged in the slow accretion phase and it 
evolves violently in the fast accretion phase. 

In the present paper, we extend and elaborate our pre- 
vious results from Paper I and analyze additional phys- 
ical quantities such as the shape, angular momentum, 
mass, kinetic and potential energies of halos. Our re- 
sults are developed in the context of the Open Cold Dark 
Matter (OCDM) scenario. However, our main conclu- 
sions are independent of the exact cosmological model 
under consideration. Since we are interested in the rela- 
tive differences between the various simulated halos, we 
shall assume here that the NFW density profile is a good 
enough approximation to a simulated halo density pro- 
file. 

The present paper is structured as following. The nu- 
merical experiments are described in § The mass ac- 
cretion history of the primary halo is presented in § EI 
and our analysis of the resulting halos is described in 
§ 14.31 Additional analysis of the structure of the DM 
halos is presented in § 0] and a general discussion in § [5] 

2. NUMERICAL EXPERIMENTS 

The common approach to study the evolution of single 
halos with high numerical resolution and from a cosmo- 
logical point of view consists of two steps. First, halos 
of interest are identified in a large cosmological iV-body 
simulation. Second, the particles which make up these 
halos are traced back to the initial conditions. The re- 
gion enclosing these particles is then re-simulated with 
a higher resolution. An alternative approach is to find 
a way of producing the requested structures in agree- 
ment with the cosmological model imposed. This can 
be done by setting up the initial density fields by CRs 
of Gaussian random field s, following the prescription of 
Hoff man fc Ribakl l(T99l|) . 

The use of CRs to set up the initial conditions as an 
input for the cosmological simulations has been applied 
to study several aspects of a structure formation. The 
CRs are also instrumental in using the data to set up 
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Fig. 1. — Model configurations. Schematic representation of our five models and their corresponding constraints. The top right labels 
indicate the mass of the constraint. The color of the labels indicate the constraint level within each circle. 



the cosmological simulations. This has been done for 
studies related to the matter distribution in the nearby 
universe fe.o...lBistolas fc Hoffmanlll998tlKravtsov et all 
120021 lRomano-Diazll2004lh However, none of these stud- 
ies have focused on the evolution of single galactic halos. 
Below we describe the main characteristics of our CRs, 
while a mathematical description of the CR formalism is 
exposed in Appendix A. 

2.0.1. The models 

We have designed a set of five different models, i.e., ex- 
periments, to probe different merging histories of a few 
times 10 12 /i _1 M o halo in an OCDM cosmology. The 
halo is constrained to have various substructure on dif- 
ferent mass scales and locations, designed to collapse at 
different times. The spherical top-hat model is used here 
to set the numerical value of the constraints. The model 
provides the collapse time of the substructures as a func- 
tion of the initial density. This is used only as a general 
and rough guide because the substructures are neither 
spherical nor isolated and hence their collapsing time 
might vary. Furthermore, the constraints used here do 
not control the experiments fully. The nonlinear dynam- 
ics can, in principle, affect the evolution in a way not 
fully anticipated from the initial conditions. Even more 
important is the role of the random component in the 
CRs. Thus depending on the nature of the constraints 
and the power spectrum assumed, the random compo- 
nent can form other significant substructures at differ- 
ent locations and mass scales. This can be handled by 
adding more constraints and varying their numerical val- 
ues. The price to pay is that the modeled entity will be 
more "synthetic." 

Table ^ presents the main characteristics of our five 
models. The first column indicates the label of the 



TABLE 1 



Model 


10 12 


Constraints (h 1 Mq) 
5 X 10 11 2.5 X 10 11 10 11 


-Zcoll 

10 10 


A 


1 






2.1 


13 


1 


2 




.3.7 


C 


1 


2 


4 


5.7 


D 


1 




6 


7.0 


E 


1 


1 


1 1 


1 8.9 



Note. — Characteristics of the CR models: the first col- 
umn indicates the model labels. The second, third — sixth 
columns indicate the level of the mass constraints expressed 
in h~ 1 Mpc. The given numbers represent the number of con- 
straints imposed at each level for each model. The last column 
indicates the collapse redshift (z co ii) of the smallest (in mass) 
constraint imposed on each model. 



model. The second-to-sixth columns show the masses of 
the different constraints (expressed in /i _1 M Q ). The last 
column states the collapse redshift of the last/smallest 
constraint in each configuration. All five models are 
embedded in a region corresponding to a mass of ~ 
W 3 h- 1 M Q in which the over-density is constrained to 
be zero — a region of an unperturbed Friedmann uni- 
verse. This spatial scale is larger by about a factor of 
three (in mass) than the size of the computational sphere. 
Therefore, this constraint cannot be exactly fulfilled, yet 
it constrains the large scale modes of the realizations to 
obey it. Model A — our benchmark model, is based only 
on one constraint of 10 12 /i _1 M©. All other models have 
substructures imposed onto this constraint as illustrated 
in Figure Q Each imposed constraint level is indicated 
by a different color and by its mass at the top right part 
of model. The position of the constraints is schematic 
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and only serves to get a visual impression of the models. 
The small-scale constraints imposed over our benchmark 
model are aimed to modify the merging history of the 
benchmark model. They change the collapse time as 
well as the number of major mergers that each model 
will pass through. In the hypothetical case, when the 
random part of the CR method does not contribute to 
the major merger history of the models, model B will 
experience only one major merger, while model C will 
have two, and model D can have more than three major 
mergers. Model E is aimed to mimic a radial collapse by 
imposing concentric (nested) constraints. 

All models have been constructed with the same ran- 
dom seed. This guarantees that the linear external 
matter distribution and the random component will 
be the same for all models Eoffrnan fc Ribald 119911 
Ivan de Weveaert fe Bertschingerlll996() . All the density 
constraints constitute (2.5 — 3.5)<7 perturbations, where 
a 2 is the variance of the appropriately smoothed field. 
The constraints were imposed on a box size of 4 /i -1 Mpc. 
In order to avoid the edge effects with the smoothing 
procedure, the constrained box has been immersed in a 
8 ft. _1 Mpc box on a cubic grid with 256 grid cells per di- 
mension. The constrained inner box was carved out and 
resulted in a box size of 4 /i _1 Mpc and 128 grid cells per 
dimension. 

2.1. N-body simulations 

Once the constrained initial conditions have been 
generated, we h ave applied the Zel'dovich formalism 
(|Zerdovicblll970f) in order to obtain the initial positions 
and velocity displacements at redshift z = 120, using 
physical rather than comoving coordinates. A sphere of 
2 ft. _1 Mpc (comoving) radius was carved out from the 
Zel'dovich evolved fields. All our constraints are totally 
embedded within this maximum radius sphere. The to- 
tal number of particles selected within this volume is 
~ 1.3 x 10 6 , resulting in a mass resolution of 3.6 x 10 6 M Q . 
The gravitational softening is 500 pc. 

We have evolved each linear density field into the 
non- linear regime, until z = (present time), us- 
ing the updated FTM-4 . 4 version of our hybrid code 
(Keller fe Shlosmanlll994t lHellerlfT995^ with the falcON 
routine IfPeTmerJ I2002T . The code is about ten times 
faster than the optimally coded iBarnes fc Hutl l|1986f ) 
tree code. The equations of motion in the FTM-4.4 
include the term describing the expansion of the uni- 
verse. The code has been suc cessfully tested ag ainst the 
Santa Barbara Cluster model ijFrenk et al.ll999|i (see Ap- 
pendix B). 

Because we are interested in following the merging his- 
tory of our models as close as possible, we have sampled 
the dynamical evolution of the systems with 165 time 
outputs spaced linearly with the expansion parameter a. 
Each halo is resolved at z = with around 1.1 x 10 6 
particles within the virial radius. 

3. MASS ACCRETION HISTORIES 

Halos grow according to a hierarchical formation sce- 
nario, in which they are assembled via merging of dif- 
ferent mass and size halos, together with a more gentle 
accretion. Several studies have found that the mass ac- 
cretion history of halos (MAH) may affect the various 
halo parameters, including R s , while maintaining their 



NFW density profiles llNayarro et al.ll997t|Bullock et alJ 
l200Tbl IZhaoet al] 12003 ITasitsiomi et all 12004 W021. 
Therefore, it is important to analyze the MAHs of our 
five models. 

W02 proposed to fit the MAHs of halos by an expo- 
nential function in the form: 

M(a) = M e- a ^/ a -V (2) 

where a is the expansion factor a nd Mq is the final virial 
halo mass. Ivan den Boschl (|2002t ) arrived to a very sim- 
ilar expression with a two-parameter f unction using the 
extend ed Press-Schechter formalism. ITasitsiomi et alJ 
l)2004|) found that in many cases Eq.|5] represents a poor 
fit to the individual MAHs halos, in particular when ha- 
los experience an intense and v iolent activity, up to the 
present time. ITasitsiomi et all proposed a more general 
MAH fit of the form: 

M(a) = a p M e^d/a-i) (3) 

where p = reduces to the W02 model. Despite the 
improvement, such a smooth fitting function cannot re- 
produce all the features in the halo evolution, especially 
the major mergers. 

3.1. Halo identification 

Groups and subgroups identification is not a trivial 
task. Several methods have been proposed which use 
different "boundin g" criteria, e.g. . Friends of Friends 
llDavis et all 119 851). DENMAX llGelb fc Bertschingerl 
I1994D. HOP fcisenstein fc Hutl li9MT SUBFIND 
ijSpringel et all 12001^ . etc. Different methods are used 
in order to isolate the halos and subhalos within the 
halos. We are interested in following the evolution of 
the main halo through the major branch along its merg- 
ing tree (see § I4.2|l . For the se purposes we use the 
publically available code HOP 5 l)Eisenstein fc Hudll998j) 
which computes groups based on the isodensity criteria. 
Next, we identify the halo centers with the densest par- 
ticle within the particular halo. Once the centers are 
found, we grow the halos by adding thin spherical shells 
and computing the total density until it reaches the value 
of A c times the mean density of the universe. Therefore, 
a halo within the virial radius will be specified as a func- 
tion of redshift as: 

Air 

M vir (z) = y A c (z) p b 0) R vil (z) , (4) 

where A c is the critical density and it is computed using 
the top-hat model, pb is the universe density and i? v ir is 
the virial radius of the halo. 

4. RESULTS 

4.1. Model evolution 

A halo of ~ 1O 12 M has formed by z = in all of 
our models, except in model B. At the present time, this 
model harbors two halos of ~ 5 x 1O 11 M which are on 
their way to merge. The final collapse time (i.e., the 
time of the last major merger) of the halos is as fol- 
lows: B, E, A, D and C, from the earliest to the latest. 
Defining the halo formation as the time when half of the 
mass is in place, results in the same order as above. In 

5 available at http:/ /cmb.as.arizona.cdu~ciscnstc/hop/hop.htinl 
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comparison, the top-hat model predicts a different order 
of collapse times. This is expected because the random 
component of the CR method can introduce other struc- 
tures of the same level as the constraints. As a result, 
the merging history and the collapse time imposed by 
these constraints will be modified (Sect ion l2. 0.11 see also 
Ivan de Weveaert fc Bertschingerlll996|) . Since we have 
warranted that our five models have the same random 
component, all of them will be affected in the same way. 

The evolution of the models can be observed in their 
density projections (Fig. EJ and the respective merging 
trees (Fig. OH- The effect of different constraints imposed 
on the halos can be noticed from their evolution on these 
diagrams. Figure El is a schematic representation of the 
halo growth and the merging history of the models. The 
size of the circles is proportional to the mass (left column 
numbers) of the depicted halos (see Section 14. 2|) . This 
representation underscores the effect of the CRs on the 
evolution of these models. The model A evolves through 
three major mergers. Since model B harbors two very 
similar halos, we only present the evolution of one of 
them (the most massive). It has two nearly simultane- 
ous major mergers of an exceptional strength, in its early 
phase. Because the halo is out of the dynamical equilib- 
rium during this time, we consider this double merger 
event as a single one. Model C frame shows the evolu- 
tion of the four constraints, their collapse into two major 
entities and formation of one single halo. In the model D, 
the six constraints and their development into two major 
halos and the final collapse can be clearly followed from 
the respective tree. The early formation of the model E 
through two recognizable major mergers is neatly shown 
in its merger tree panel. 

Figure El exhibits density maps in the co-moving co- 
ordinates of the models for five different epochs taken 
from the respective merger trees. The first column refers 
to model A, the second — to B, the third — to C, the 
fourth to D, and the last one — to model E. The 
color scale is proportional to the logarithm of the density 
field. The first row shows the different initial conditions 
between the five models. The last row shows the final 
state of the halo. Shown are the random frames after 
the last major merger that each model experienced, and 
not the last snapshot of our simulations. The intermedi- 
ate frames illustrate different dynamical epochs for the 
halos. In model A, the frames depict the moments prior 
to two major mergers. Model B frames clearly show the 
two imposed constraints — the system, although grav- 
itationally bound, does not merge at the present time, 
but it will do this in the future. Models C and D show 
the early collapse of the small constraints, the forma- 
tion of the two major clumps and their way to collapse 
into a single halo. Model E (the first model to collapse 
into a ~ 10 12 M Q halo) shows a more radial filamentary 
structure around it, a direct consequence of the imposed 
constraints. The amount of substructure changes from 
model to model, with the model D exhibiting the richest 
substructure. 

4.2. Merger trees 

Since we have identified all halos for basically each of 
the output times of our simulations (we have found halos 
for all models from z ss 20 onwards), we can construct 
detailed merger trees along the path of the main pro- 



genitor. Using the main halo at a given time, we trace 
back the identified particles to the previous time step and 
correlate them with the most massive halos at this time. 
The progenitor is identified with that halo which contains 
at least 50% of the particles of the parent halo. We also 
follow the trajectories of their respective centers. This 
will assure that we follow the right predecessor. Once a 
branch has been chosen along the halo's tree, we follow 
it up until this branch opens again and repeat the same 
procedure. Note that in this way the followed halo is 
not necessarily the most massive halo at all times in the 
simulation, but rather most of the time. 

Figure El shows the merger trees for our five models. 
The size of the circles is proportional to the logarithm 
of their respective masses (indicated by the numbers lo- 
cated by the left columns of each panel) at the respective 
times (indicated by the right columns labels). The de- 
gree of accuracy of the trees depends upon the number 
of identified structures at each time step. In general, our 
merger trees consider all identified structures per time 
step, but for the sake of clarity we have only plotted the 
relevant halos to illustrate the evolution of the models. 

Once the merger trees have been elaborated, one can 
construct the MAHs straightforward. The red dots 
in Figure El represent the measured masses normalized 
with respect to their respective final halo masses. The 
blue solid lines represent the best fit according to the 
W02 model, while th e dashed black lines — to the 
iTasitsiomi et alJ l)2004ft model. The fits were obtained 
by x 2 minimization. In general, the ITasitsiomi et all 
model represents a better fit to all models than the W02 
one. However, both methods fail in reproducing the early 
parts of the MAHs. The deficiencies of the W02 model 
are emphasized in models C and D which have a more 
violent his tory than the res t of the models. On the other 
hand, the ITasitsiomi et al.l fits do not represent the full 
MAH of such halos, see for example the large deviations 
at a ~ 0.6 for both models (the time between two major 
mergers) . It is clear that neither fits reproduces the MAH 
of a halo throughout its entire history. At best such fits 
can accurately describe the MAH within one single pas- 
sive evolution phase. Having stated that, we have used 
Model B, which has the longest quiescent phase, to test 
the quality of the analytical fits. 

4.3. NFW analysis 

For all halo models, we fit the NFW profiles as a func- 
tion of time and thus follow the cosmological evolution 
of their parameters (R s and p s ). Many factors may af- 
fect the fits of analytical profiles to the observed (sim- 
ulated) ones: the choice of binning, the merit function, 
the weights assigne d to the data points, the range of radii 
used in the fitting ijTasitsiomi et alJl2004(l . We have di- 
vided each halo into spherical shells equally spaced log- 
arithmically until 0.6i? v i r , in order to give more statisti- 
cal weight to the inner regions of the halos. We have 
chosen this radius since fitted profiles do not change 
substantially beyond this radius. In case that a mas- 
sive subhalo is nearby (we define such halos as those 
with mass > 0.3M v j r ) the density profiles are non- 
monotonic and bumps at the data point distributions are 
present. We avoid such bumps by performing fits until 
mm(0.6i? v ir 7 0.5dh), where c?h is the distance to the given 
subhalo. We have set the minimum fitting radius equal to 
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Fig. 2. — Time evolution snapshots for the five models. Models are arranged from left to right columns: A, B, C, D and E. Each frame 
represents the logarithmic density field weighted by the mean density at the given time. The top left labels in each frame indicate the 
rcdshift of a given snapshot. 



2e (where e is the gravitational softening imposed in the 
JV-body force calculation) in order to avoid two-body in- 
teractions. This minimum radius is well within the range 
where a density profile for the characteristi c of our mod- 
els can be considered to be "resolved" ( e.q., lBinnevl2004t 
IDiemand et al.l2004tlReed et al.l2005a|) . The fitting pro- 
cedure is performed by weighted x 2 , where the residual 
in a given shell is normalized by its own density. 

The upper panel of Figure shows the NFW fits for 
12 different time outputs (indicated by different colors 
and symbols) of model C. The bottom panel presents 
the same fits but normalized by R s and p s respectively. 
Symbols indicate the data measurements while the con- 
tinuous lines the best fitted models. Most of the fits are 
in good agreement with the measured data. Such fits cor- 
respond to the steady/quiet epochs in the halo evolution. 
The fits that do not follow the observed data correspond 
to those epochs where the halo is passing through a vio- 



lent episode (major merger). 

4.4. Virial quantities 

By following the evolution of the halo virial quantities 
(i? v ir, Afvir), we arc able to distinguish, as a function of 
time, the main large scale characteristics between the dif- 
ferent models. Figure [S] shows the evolution of the virial 
radii for the different models (solid lines). In the bottom 
right-most panel, all profiles have been superimposed for 
a better comparison. Note that the virial radius grows 
almost linearly with a during the quiescent phases in all 
models. The sudden increases correspond to the epochs 
where violent activity takes place indicated by the square 
brackets at the bottom of each panel. The width of the 
brackets shows the duration of the violent event. Despite 
the fact that the various models have different violent 
epochs at different times, all of them (apart from model 
B) seem to converge to the same virial values once their 
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Fig. 3. — Schematic merger trees for our five models, 
of the halos at the given times. 



The right columns indicate the redshift frame. Circles are proportional to the mass 



last major merger has taken place. This happens since 
the larger constraint has been set to form a ~ 10 12 M Q 
halo. The violent activity is more noticeable in the virial 
mass (Figure continuous lines) where the large jumps 
indicate a substantial mass component addition to the 
main progenitor. During the quiescent epochs, mass is 
only deposited very slowly through accretion and minor 
mergers that are being tidally disrupted outside R s . This 
can be observed as a very gentle increase in the mass tra- 
jectories between the major mergers. 

4.5. Inner characteristic quantities: R s , M s , c, etc. 

Figure El also exhibits the evolution of the inner ra- 
dius R s (broken lines), and the concentration parame- 
ter c (dotted lines). Subject to a ~ 10% jitter, R s is a 
growing step function. In other words, R s grows only 
at the moments when major events take place (indicated 
by the bottom square brackets) and it remains constant 
along the quiescent phases. The horizontal lines which 
overlie the distributions are the average of R s between 
the major events. Note that the R s distributions follow 
closely the mean lines for the whole ranges. The ampli- 
tudes of the jumps are proportional to the degree of "vio- 
lence" of the encounters, in other words, to the change in 
the kinetic energy that the system experiences (Paper I 
and Sec. EH). The large jumps and discontinuities in R s 
are a clear indication of major mergers taking place. In 
the bottom right-most panel, all R s trajectories can be 
compared more closely. All models reach approximately 
the same value after their last major merger (apart from 
model B), independently of the number of major merg- 
ers and the epoch they took place. It is interesting to see 
that model B has the largest R s among all the models. 



The corresponding concentration parameters c (Fig- 
ure [SJ) display an erratic behavior at early epochs, e.g., 
a ~ 0.1. However, starting with a ~ 0.2, the trajec- 
tories become more coherent. Different model R s and 
i? v ir exhibit the tendency to converge after their last ma- 
jor merger, while c does not (the bottom right panel of 
Fig.lHJl. The jitter present in c is mainly dictated by the 
R s behavior. It has been claimed that c grows linearly 
with time (W02). However, since we resolve the step- like 
evolution of R s and the discontinuities in i? v ir, c appears 
to grow linearly only during the quiescent phases. We 
show this clearly in Fig. |H| where the black continuous 
lines represent the best linear fits for the given quiescent 
phases. The numbers shown within each frame represent 
the best fit parameters for the slopes, m, and zero-points, 
s. Note that the slopes decrease after each major event. 
The intermediate slopes for Models A and E are unreli- 
able and were not fitted because the halos do not have 
enough time to relax between these major events. 

Model B has been assembled earlier compared to other 
models, but it possesses the smallest c at the present 
time. This behavior appears to be in contrast with 
the previously reported results fi.e. lBullock et al"ll2001bL 
W02) which claimed that the halos assembled early, when 
the universe has a larger density, are more concentrated 
and have a larger c parameter. However, R s and, there- 
fore, c are related not only to the density of the universe 
at that time, but also to the intensity of an event dur- 
ing the formation epoch of this particular halo (Fig. |H1 
see also Paper I). Since model B has experienced the 
strongest event in terms of the deposited kinetic energy 
K s within R s (see Fig. lilt , it has the larger R s and, 
therefore, the smallest c. 
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It is interesting that both models C and D seem to 
depart from the general trend observed by the rest of 
the models at a ~ 0.7. This corresponds to the epoch in 
which both models went through their last major merger. 
After these events took place, the c trajectories lower 
their amplitude and join the general trend. The growth 
of c seems to be nearly linear with respect to a during 
the quiescent phase, indicating that its behavior is dom- 
inated mainly by i? v ir (see also Fig- EJ • 

The mass M s is computed by counting particles within 
this radius and hence it is a top-hat mass. The M s be- 
havior is presented in Figure [7| (dashed lines) for all five 
models, and is very similar to that of R s — it remains 
constant during the quiet phases and grows only during 
the violent events. The horizontal lines indicate (as in 
the case of R s ) the mean M s within the given time in- 
tervals. Similarly to Fig. El the bottom brackets in each 
frame indicate the time and duration of the major merg- 
ers. The model B also has the largest M s , as noted at 
the bottom right-most panel where all M s trajectories 
are depicted. At the time of the violent events associ- 
ated with the core growth, the system is not in a virial 
equilibrium and so the density profile is undefined caus- 



ing the erratic appearance of spikes in R s and M s (e.g., 
Figs. EH3 and © . 

4.5.1. General trends in R s 

The R s behavior detected in Paper I and in the pre- 
vious section has been only addressed from the 'typical' 
NFW fitting procedures. As has been pointed out in Sec- 
tion ^ there is a controversy whether the inner slope is 
— 1, steeper or shallower, or it is not universal at all. This 
raises the question whether the R s behavior is driven by 
the change in the inner slope (a mere effect of the NFW 
fitting procedure) or by a real change in R s . To answer 
this question, we fit the halos using three alternativ e den- 
sity p rofiles, i.e., by applying the generalized law ijZhaol 
Il99ffl : 

^ (r/R^il + ir/R^)^' 

where (a, 3, 7) = ( 1, 3, 1) for a NF W, (1.5, 3, 1.5 ) for a 
IMoore et all lll999l) . (1,3,1.5) for a l.Ting fc Sutol ll200"l 

and (a, 3, 7) for a more general density profile. In gen- 
eral, the fitting procedure varying the parameters a, (3 
and 7 leads to degeneracies between the parameters (see 
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iKlvpin et all2001tlTasitsiomi et alJl2004fl . With this ap- 
proach, it is possible to find "excellent" fits using param- 
eter values far from what have been found for that kind 
of halos. Therefore, in order to get physically meaningful 
fit estimates one has to constraint the parameter range. 

Figure [§] shows the evolution of i? s 's obtained from 
four different fitting profiles indicated by different colors. 
Note that all profiles show the same behavior, R s remains 
constant during the quiet phases and increases abruptly 
during the violent epochs. It is expected that the applica- 
tion of different density profiles will resu lt in differing val- 
ues o f i? s , which can be correlated fe.p.. lTasitsiomi et al.l 
|2004|) . The general fits (red lines) are very close to the 
NFW fit (black lines). This effect could be due to our 
u se of NFW values t o constrain the general fitting. 

iReed et"aTI l)2005a|) found that the inner slope does not 
decrease in time, but rather stays constant during the 
quiescent phase, by applying a semi-general fitting to the 
halo evolution. We can apply this observation to each of 
the phases and conclude that the slope also changes only 
during the violent events. Combining these two effects, 
one expects that the NFW R a displays a similar behav- 
ior to R s obtained from other fitting procedures; Fig. [5] 
shows clearly this point. Because we observe that R s 



grows as a step function even for the general fit, this 
indicates that not only a change in the inner slope is 
involved, but rather a real behavior/change in R s . 

4.5.2. The characteristic density p s 

The characteristic density has the opposite behavior to 
R s — it is a decreasing function of time, so their prod- 
uct, p s R s , appears nearly constant with time for all mod- 
els and across the violent and quiescent phases. As i? s 
also p s is affected by the 10% jitter. At major merg- 
ers it drops considerably and fluctuates until it reaches 
a plateau where it remains constant. At the next major 
merger it drops and reaches a new, lower level plateau. 

We have investigated the possibility of a correlation be- 
tween R s and p s within the jitter of the quiescent phases. 
Naively one expects that R s and p s in the jitter will pro- 
duce a scatter diagram. However, Figurc lTUI shows that a 
strong correlation exists between both quantities: larger 
R s correspond to smaller p s . Moreover, a one-to-one cor- 
relation exists between R s and p s across the full range of 
R s . The different colors represent the various quiescent 
episodes that the halos experience through their MAHs. 
The black symbols correspond to the violent episodes 
and to the early assembly time of the halo. The scat- 
tered points located above and below the main trends 
mark the violent events. Model B displays a single phase 
since this halo formed at an early phase and remained 
quiescent for most of its evolution. Models C and D, 
which display three well identified quiescent phases, have 
three well defined point distributions separated from each 
other. Models A and E present very similar distributions. 
We have performed linear fits to the logarithmic distri- 
butions for all models of the form 

log p s = m log R s + b . (6) 

in two fashions. First, by including all plotted points; 
second, by separating the quiescent phases and fitting 
each. The upper legends at each panel indicate the slope 
of the individual fits (m) and their zero points (b). Note 
that all the zero points lie about log p s ^ 18.2. The slope 
averaged over all models is m = — 1.37 ± 0.03. Even the 
slope constructed from different epochs (given by differ- 
ent colors in Fig. I10[l shows very little scatter around 
the average value. This value appears to be indepen- 
dent of particular halos and of their evolution. Further- 
more, the points corresponding to the later time (i.e., 
red color) are consistently located in the tail of their re- 
spective distributions. This color order is a direct conse- 
quence of R s and p s monotonic evolution discussed ear- 
lier. The measured jitter (10% or so) in both R s and 
p s can be recognized now as the individual trends (col- 
ors) for each model. This implies that a clear order 
exists in the halo variations. Specifically, we find that 
yO s ~i? s _159±015 during the jitter at early quiescent, i.e., 
a < 0.5, epochs in all models, while p s ^i? s _1 - 39±0 04 for 
later times. When translated to M s — i? s correlation, this 
amounts to M s ~i? s 141 and M s ~ii s 1 ' 61 , respectively The 
simplest interpretation of these correlations is that they 
are driven by density fluctuations which originate at radii 
between the cusp (i.e., slope = -1) and R s (i.e., slope = 
-2). 

Lastly, the shifts between different colors in Fig. ITU1 
come from violent events which deposit kinetic energy 
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Fig. 6. — Evolution of the virial radii R vlI (solid lines), characteristic radii R B (dashed lines) and concentration parameter c (dotted 
lines). Each panel refers to a different model. The bottom right- most panel compares directly the i? v i r and R B quantities. 




Fig. 7. — Evolution of the virial mass M v i r and the characteristic mass M B (in units of the total mass, Mtot) for all models. The bottom 
right-most panel presents all mass trajectories. Different colors correspond to different models. 



and mass within R s and separate the quiescent evolu- 
tion (see section 5. 5). Our resu l ts are in agreement with 
those reported by iZhao et all l)2003j) , who found that 
M s ~i? s 1 ' 44 holds in the "slow accretion phase" (gentle ac- 
cretion), and M s ~i? s 192 in the "rapid accretion phase." 
The difference results from our ability to separate the vi- 



olent and quiescent phases — we have excluded the black 
symbols in Fig. from the fitting, because they corre- 
spond either to the very early assembly or later violent 
times. This reduces the scatter in our fits below that of 
Zhao et al. 
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Fig. 8. — Linear plot of the evolution of the concentration pa- 
rameter, c, for the five models. The individual curves has been 
smoothed over ten points. The continuous lines represent the best 
linear fits for various quiescent phases. The fitting has been avoided 
for short phases because of the associated errors. The enclosed top- 
left parameters indicate the best zero-points, s, and slopes, m. 
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Fig. 10. — The characteristic density p B as a function of the 
characteristic radius R B . The colors indicate the different quiescent 
phases that each halo passes through. The individual slopes, m, 
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for every model (see text for further details). The fitting has been 
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4.6. Kinetic and Potential energy and the virial ratio 

We have computed the internal kinetic energy K of the 
main halo along its merger tree. For each time output, 
we have computed the kinetic energy, K r = ^T=i \ m ^ G \i 
where the summation goes for all objects enclosed within 
a radius r, n r , and o~\ represents the velocity dispersion 
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Fig. 11. — Kinetic energy for all models as a function of time. 
The blue distributions represent the characteristic kinetic energy, 
while the red ones the virial distributions. Labels are the same as 
in the previous figures. 



of a given object i of mass mj. 

Figure ^] shows the evolution of the internal kinetic 
energy computed within the radii R s (blue lines) and 
JJvir (red lines). As in the case of R s , K s also increases 
as a step function, i.e., only at the major events and 
it remains constant during the quiescent phases. On the 
other hand, i^vir grows in two different ways, it grows 
suddenly at the violent phases as K s does, while during 
the quiescent phases it grows in a very gentle way. 
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From Fig. ^] one can realize that K s of model B has 
the largest value of all models while its K viY is the small- 
est from all models. This is a mere reflection of their mass 
accretion histories. Furthermore, the growth of model B 
energy curve is the most abrupt of all models, occurring 
at an early redshift [z ~ 9). This halo went through a 
very violent and rapid formation epoch at this early z. 
This violent formation event resulted in the largest i? s 
in our sample, although it is the less massive of the five 
models. 

The potential well of the halos is also co nstructed from 
their respective MAHs. iZhao et"aTI l^OOffi found in their 
simulations that the potential well is mainly build up 
in the so-called fast accretion phase, while remaining al- 
most constant during the slow accretion phase. Figure lT^l 
shows the evolution of the potential energy (0) for our set 
of models. This <f> is computed by including all particles 
within a given radius, namely R s and i? V ir- Since there 
is no mass accretion within R s in the quiescent phases, it 
is expected, therefore, that the potential </> s remains con- 
stant during such phases and only deepens at the major 
events when mass is accreted within R s . In the case of 
the potential within i? v ir, where the halo growth is due 
to violent and gentle mass accretion, </> v ir uncovers these 
two phases as well. It sinks substantially at the major 
events and it deepens in a very gentle way during the qui- 
escent phases. This general behavior indicates that the 
potential wells are only formed and substantially modi- 
fied at the violent events. The main difference with the 
IZhao et "ail l|20f).1fl analysis is that they only recognize a 
single early fast accretion period, while in our models the 
later quiescent phase is intermitted by the violent events. 

The way energy is acquired by the halos and the de- 
tails of the formation of the respective potential wells 
reinforces our view of the halo growth process. Energy is 
deposited into the halo via gentle and violent mass accre- 
tion. The latter affects the whole halo since both R s and 
i? v ir grow during such events, while the former affects 
only the exterior part of the halo. Gentle accretion and 



minor mergers deposit their mass and energy within this 
outer region. Energetic events such as those produced by 
the major mergers are able to reach the halo core, since 
they are more strongly gravitationally bound. Part of 
their mass and energy are deposited inside the halo core. 

The interplay between the kinetic and potential en- 
ergies within a given halo can be better understood in 
terms of the virial ratio defined as Q r = —2K T /<f> T , where 
the subindex r denotes the radius at which this quan- 
tity is computed. At the virial radius, Q r should be 
close to unity for systems in virial equilibrium and nearly 
isolated. Because virialized structures are bordered by 
infalling and outgoing material, the convergence of Q T 
to unity at Ryw should not be fully expected (see e.g., 
iMaccio et all l2003 LShaw et alJ l20f)l . However, since 
our models are basically constructed in isolation, after 
their respective last major merger, Q v ir^ 1- Indeed, 
this can be observed in Figure EH where the solid lines 
represent the Q v ir quantities for each model. The pro- 
nounced spikes along each track represent the time of 
the major events and the considerable departure from 
the virial equilibrium. On the contrary, Q s is not ex- 
pected to be close to unity, since this region experiences 
"pressure" from overlying halo. Note that during the 
quiescent phases Q s remains constant and, as in the case 
of R s , jumps to another stationary trajectory when the 
major events take place. 

4.7. Shapes 

We have analyzed the shapes of our halos as a func- 
tion of time by using the inertia tensor, for R s and i? V j r 
respectively, using the unweighted moment of inertia de- 
fined as /ij = J2n=i x i,nXj,n, where x- u Xj represent the i, j 
Cartesian components of the position of a particle n and 
the summation extends over all pa rticles enclosed within 
the given r adius r {= Rs, i?™)- lAllgood et all pool 
showed that the weighted and unweighted moments of 
inertia used to compute halo shapes differ only by 10%. 
Our aim here is to compute the halo shape evolution 
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Fig. 14. — Axes ratios for our five models. For each model we 
present s (= c/a, dashed lines) and q (= b/a, solid lines) ratios for 
R B (blue lines) and R v i r (red lines). 
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at various radii as a function of their merging history. 
Therefore, we divide the halos into spherical shells. An 
alternative approach is to compu te the shapes within iso- 
surfa ces of density or potential ijBerentzen fc Shlosmanl 
|2006|) . Since so far we used quantities defined within the 
R s and i? V ir spherical distributions, we choose the former 
approach for consistency. 

We compute the eigenvectors and eigenvalues of the 
moment of inertia matrix and define the axes of the el- 
lipsoid as a, 6, c = -\/(Aa, At,, A c ), where Ai corresponds to 
the ith eigenvalue, and they are arranged in the descen- 
dent order (i.e., A a > Ab > A c ). These axes are nor- 
mally described in terms of the ratios s = c/a, q = b/a 
and p = c/b. 

Figure ^] displays the shape evolution of R s and i? v ir 
indicated by the red and blue colors respectively. The 
dashed lines correspond to the s parameter of both radii, 
while the solid lines — to the q number. Clearly, all mod- 
els are prolate within R x . This result is in agreement with 
th e analysis of | A llgood e t alJ 12006|) at small radii, and 
of IBerentzen fc Shlosmanl l|2006j) . Within R s the halo 
shapes are unperturbed and only change at the major 
mergers. On the other hand, the i? v i r shapes change dur- 
ing the quiescent epochs as well, mainly because gentle 
accretion is not isotropic. The s and q at i? v ir indicate 
that there is a tendency for the halos to become less tri- 
axial with time. 

4.8. Angular momentum considerations 

One of the drawbacks of the present simulations is the 
lack of the external tidal field presence in our simula- 
tions. Although our simulations have accounted for the 
tidal field within the original box of 8 /i _1 Mpc, this vol- 
ume does not suffice to include the large scale tidal con- 
tributions. Nevertheless, this affects only the final phases 
of the formation of the halo. At earlier times, the main 



halo is surrounded by other halos of about the same mass 
that torque it. Therefore, the measured angular momen- 
tum will be mainly produced by the direct interactions 
of the halo with those substructures created with the 
constrained random field. Any differences between the 
models will be attributed to their different merging his- 
tories. 

The angular momentum can be expres sed in term o f 
the dimensionless spin parameter (e.g., iPeeblesI Il969j) 
A = J\E\^ 2 /(GM 5 ^ 2 ), where J, E and M are the 
total angular momentum, energy and mass of the 
system, and G is the gravitational constant. The 
value of the spin parameter roughly corresponds to 
the ratio of the angular momentum of a halo to that 
neede d for rotational support (see e.g., iPadmanabhanl 
119931) . Typical values of the spin parameter of in- 
dividual halos fro m A-body simulations are in the 
range [0.02. 0.111 jBames fc Efstathioul Il987t iRvdenl 
1 98811 Warren et alll992tlSteinmetz fc Bartelmaimll995l 
Cole fc Lacevll199aiCard n erl200lUBullock et 3,1.112(101 al). 
The spin parameters in A-body simulations are inde- 
pend e nt of the cosmological model iBarnes fc EfstathiovJ 
19871 iWarren et all Il992l iLemson fc Kauffmannl 1199a 
Cardnerl 120011). halo env ironment or halo mass 
ijLemson fc KauffmannlH 999(1 . The only correlation that 
has been found is with the time of the last major merger: 
halos tha t had a recent m ajor merger have larger spin pa- 
r ameter (lCardnerll2001l) . 

iBullock et alJ l|2001a|) proposed an alternative and 
more practical way to compute the spin parameter A' = 
J/(\^2 M V R), where J is the angular momentum within 
a spherical volume of radius R with mass M and V is the 
halo circular velocity at radius R, V 2 = GM/R. This 
definition reduces to the standard one when measured at 
the virial radius of a truncated singular isothermal halo 
(|Bullock et al.ll2001a|) . In A' definition, J is defined as 
J = J27=i r i x TO i v i- We have adopted the IBullock et all 
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definition, hereafter A. Figure ^] shows the evolution 
of the spin parameter computed at the two usual radii, 
i? V ir and i? s , A s and A V i r (blue and red lines respectively). 
The sudden changes in A v i r are clearly correlated with 
the major events involving each halo. A v j r exhibits the 
tendency to decline d uring the quiescent epochs (see also 
IVitvitska et"afll2002j) . This is a consequence of fast in- 
crease in i? v ir~ i during these epochs, which makes the 
denominator in A to grow faster than J for any reason- 
able scenarios of the halo growth. Models A and E have 
the lowest virial spin parameter because most of the ac- 
creting material has already joined them. Model B has 
the largest A V i r by the end of the simulations — it has the 
smallest main halo and a substantial amount of material 
is locked in its companion and around it. 

The behavior of A s differs somewhat from that of A v i r - 
It confirms that the halo core grows only via accretion 
of major mergers and minor mergers do not play a role. 
Under these conditions, one expects that A will be con- 
served. The blue lines of Fig.lT^lshows precisely this type 
of behavior. The spin parameter decreases only at the 
major events and remains constant during the quiescent 
phases in most of the models which exhibit the same A s 
value. However, Model B has the lowest A s , which is due 
to a smaller number of major mergers it has experienced. 

This analysis shows that although the models miss 
large scale torques, this does not represent a consider- 
able problem for their evolution. While the models have 
been constrained to the same initial sizes and masses, 
no constraint has been imposed on the spin parameter. 
Because the number of major mergers is the same in all 
models except the Model B, the final spin parameter dif- 
fers little between them. 

5. DISCUSSION AND CONCLUSIONS 

We have investigated the effect of a different assembly 
history on the final structure of the DM halos. We have 
employed the Constrained Realization method of random 
Gaussian fields in order to create the initial conditions 
for a 10 12 h~ 1 M Q halo, which has been evolved subse- 
quently by means of the iV-body simulations. Five dif- 
ferent variants of the same final halo have been simulated 
and analyzed using a sufficient mass resolution to iden- 
tify the halos and subhalos at early redshifts (z > 10), 
and a sufficient time sampling in order to closely follow 
their dynamical evolution. The model evolution has been 
quantified in terms of parameters which characterize the 
NFW as well as other, non-NFW density distributions. 

The evolution of a given halo can be characterized by a 
number of quiescent phases of a slow evolution intermit- 
ted by violent episodes of major mergers. We find that 
the inner halo is in a state of a dynamical equilibrium 
during the former phases, and its density is well approxi- 
mated by an NFW profile. Furthermore, the NFW char- 
acteristic radius R s appears to be the best gauge of the 
dynamical evolution of a halo. It remains constant in the 
quiescent phases and it grows discontinuously during the 
violent episodes. Other variables that characterize the 
inner structure, i.e., within i? s , behave in a similar way 
in their growth or decline, e.g., p s , M s , A s , etc. Between 
the major mergers, the halos evolve very passively. This 
means that a gentle accretion does not influence the be- 
havior of the characteristic quantities. The minor merg- 
ers only contribute to the gentle growth of the external 



halo and to a decrease of the spin parameter. We note 
that the product p s R s is nearly constant in time showing 
only a very slight decline. The significance of this will be 
discussed elsewhere. 

On the other hand, the virial parameters, e.g., i? v i r , 
M v i r , A v ir, etc., exhibit a monotonic growth or decline 
during the quiescent phases, in ac cordance with the sim - 
ple analytical relation (e.#. ; W02, iBullock et al"]l2001b|) . 
During this time, the virial radius i? v ir shows a nearly 
linear growth with the expansion parameter a, and the 
mass accretion history, namely M v \ r , is fairly well de- 
scri bed by the fitting formu la of IWechsler et all l|2002f) 
and iTasitsiomi et all |2Q04) which exhibits an overall 
slowdown for the later times. However, these simple re- 
lations hold in the quiescent phases only, and, in general, 
cannot be extrapolate from one phase to the other. Dur- 
ing the violent episodes, the virial parameters change 
discontinuously. Therefore, one cannot fit the entire evo- 
lution of a halo, consisting of a a series of violent and qui- 
escent phases, with a single analytical expression. The 
occurrence of a few generations of violent phases is typ- 
ical for the CDM cosmogonies rather than an anomaly 
and any theory which aims at modeling the evolution of 
halos should account for that. 

The above discussion implies that the concentration 
parameter c closely mimics the behavior of i? v i r during 
each quiescent phase by showing a nearly linear growth 
with a. This growth, however, slows down with each 
consecutive quiescent phase. The value of c at the on- 
set of each slow phase appears to depend both on the 
structural change in the halo that underwent the violent 
merger (i.e., change in R s ) and on the degree of violence 
of that merger. 

We note that the trends observed here do not pertain 
only to a NFW fit but remain when various fitting formu- 
lae for the density profiles have been applied. The actual 
numerical values of the structural parameters (i.e., R s , 
R vlr , etc.) have been found to depend on the exact fitting 
procedure, but the general trends did not change. There- 
fore, our conclusions reflect the robust characteristics of 
the halo evolution. 

The original goal of this work was to impose different 
constraints on the mass distribution within the computa- 
tion box with the same total mass. Whether this would 
lead to a diverging evolution was one of the main objec- 
tives of the current project. One of the aspects was to 
obtain a different number of major mergers within the 
box. The crucial question was to what extent the final 
halo properties depends upon its assembly history. Our 
results have shown that the imposed constraints on the 
mass distribution had only a partial effect on the initial 
conditions because of the random component. As a re- 
sult, the number of major mergers after z ~ 10 was the 
same (i.e., three mergers) in all the models, except in 
model B (two overlapping mergers or one). Hence the 
initial conditions did not possess a sufficient scatter to 
form the final halos characterized with wide range of in- 
ternal structure. Nevertheless, clear trends emerge. 

First, one is unable to reconstruct the evolution of the 
halo without accounting for the amplitudes of the discon- 
tinuities triggered by the major mergers. Their overall 
effect is non-additive, in the sense that larger fractional 
energy inputs, AK S /K S , during the mergers result in the 
non-linear growth of R s and other structural parameters 
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as well. 

In order to quantify the violent phases of halo evo- 
lution, we introduce the 'violence' parameter, Q V i = 
AK S /K S , defined in terms of the fractional kinetic en- 
ergy deposited within R s in a major merger. Assuming 
virial equilibrium this parameter measures the deepening 
of the gravitational potential in a violent event. We de- 
note the change in the characteristic parameters across 
the major merger event by AR S /R S = (i? s2 — i? s i)/i? s i, 
where subscripts '1' and '2' refer to 'before' and 'after' 
the event. Fig. 2 of Paper I shows that the fractional 
change of R s in terms of the fractional change of K s are 
approximately related by 

where typically < Q v k> ^ 1, but it also can be some- 
what larger than unity in principle. Thus the change in 
R s varies in a non-linear fashion with Q v io- Therefore, 
the evolution of R s depends both on the number of the 
violent phases and the magnitude of each one. This im- 
plies that the evolution of R s cannot be reconstructed 
from integral quantities but it depends on its detailed 
merging history. 

Second, we compare our models at four benchmark 
rcdshifts, namely z = 0,1,2 and 3. For this purpose, 
we only use the primary halos with the same virial pa- 
rameters, i.e., of the same virial mass. Figure shows 
that models B, C and D possess similar mass halos at 
z rsj 1 — 3. However, the structural differences between 
these halos at higher redshifts are more substantial than 
at z = 0. For example, The NFW scaling parameters 
R s exhibit much larger scatter at earlier times. They 
converge towards z = 0, except for model B. The latter 
has the largest R s despite going through an early pro- 
longer major merger — an indication that the energetics 
of mergers characterized by Q vlo can be as important to 
the halo's structural evolution as the number of merg- 
ers (Fig. EJ. At z ~ 1 models B, C and D observe the 
relation R^ ~ R® ~ 0.5i?f . This is also the case of 
models D and E at z ~ 5., for which Rf ~ 0.6Rf . These 
examples show how the merging histories of otherwise 
similar objects can affect their NFW structure, i.e., the 
R s related quantities. 

Third, the value of the concentration parameter c has 
been demonstrated to depend not only on the formation 
time of the halo but also on the degree of violence in 
its merger events, Q v i . Hence, the halos forget the ini- 
tial conditions to some extent. This trend is expected to 
become more obvious with the increase of the computa- 
tional box, and, hence with the associated increase in the 
number of major mergers that the halo goes through. 

Lastly, we comment on the effect of the environment 
on the growth of the halos. The virial parameters in our 
models grow nearly linearly (on linear plots!) during the 
successive quiescent phases. A closer look on Figs.EJand 
|H1 reveals that this growth tends to flatten at later times 
in all models except model B which stays linear since its 



early merger. This latter model is the only one where the 
largest constraint is still collapsing at z = 0. Hence halos 
which are located in rich environments are expected to 
continue growing even at present time while for those 
which are located in voids, the growth will saturate. 

The emerging picture from the evolution of character- 
istic masses and radii and of virial masses and radii con- 
firms that only the major merger s are able to penetrat e 
the core, as has been argued bv iSver fc White! l)1998|) . 
The minor structures seem to be stripped completely in 
the external halo before they could reach the inner core, 
and so the mass of the external halo region grows during 
the minor mergers and the gentle mass accretion. This 
relatively straightforward picture will be complicated if 
baryons are included. Dissipative processes associated 
with them will lead to the formation of a substantially 
more bound systems, either purely baryons or mixed with 
the DM. Such systems will be able to survive during very 
un-equal (minor) mergers and contribute to the evolution 
of the halo core. 

The main limitation of the present simulations is the 
small size of the computational box. In none of the mod- 
els the halo is embedded within the typical cosmic web 
and it evolves rather as an increasingly isolated halo. The 
size of the box limits the effect of the tidal interactions 
which induces the halo's angular momentum. Neverthe- 
less, the resulting spin parameter of the halo in all sim- 
ulations, A w 0.01, is still within the scatter measured 
in cosmological simulations. One can conjecture that to 
the extent that the large scale tidal field affects the in- 
ternal structure of a halo, the increase in the size of the 
computational box might lead to a stronger dependence 
of the halo structure on its merging history and hence on 
the imposed constraints. 

The five models considered here provide five different 
possibilities for the merging history leading to the forma- 
tion of essentially the same halo. Thus, only one main 
halo is considered here, based on one given realization 
of the primordial perturbation field. More realizations 
are needed to construct an ensemble of halos so as to 
calculate the evolution of the averaged quantities. It fol- 
lows that actual numerical values of quantities of interest, 
such as the concentration parameter, R s , and so, cannot 
be considered as representative. Yet, the different evolu- 
tionary trends pointed out here provide an accurate de- 
scription of the evolution of DM halos in the CDM-likc 
cosmologies. 
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APPENDIX 

CONSTRAINED REALIZATIONS FORMALISM 

Consider the primordial cosmological density field 6(r) on which a set of linear constraints are to be imposed. The 
constraints consist of a set of positions, mass scales and the value of the 6 field at the specified locations, appropriately 
smoothed on the mass scales, namely {r Q , M a , A a } 1 M , where r Q is the location, M a the smoothing mass scale 

and A Q is the value of the smoothed density field of the a-th constraint. The smoothed field is define by a convolution 
with a Gaussian kernel, 

A fl( r ) = / d3fc ^cxp 



k 2 R 2 



(2tt) 3 

where 5^ is the Fourier transform of the primordial perturbation field, 
Gaussian smoothing, R is related to the smoothing mass scale by 



S, and R is the linear smoothing scale. 



(Al) 
For a 



R = 0.64 



M 



-1 1/3 



(A2) 
The 



(47r/3)/9 cri t^i 

where p cri t is the critical cosmological density and f2 m is the matter density parameter (Bardecn et al. 1986) 
c onstraint on the mas s scale M is referred to as on the co-moving linear scale R. 

iBertschinsrerl l)1987j) proposed that a constrained field can be composed from the sum of a mean field (fixed by the 
values and forms of the constraints) and a residual field. This last field ad ds the random componen t to the mean field. 
A constrained realization of the unsmoothed S field has been obtained bv lHoffman fc Ribald l)1991[) : 

5 CR (r) = ~S(v) + (S(v) A Rq ) (a Rq A Rg y\A Rg - A % ). (A3) 

Here, <5(r) is a random realization of the S field and Ar q is a mock constraint obtained from the random realization, 



A R (r Q 



(2tt) 3 



d k 6^ cxp 



k 2 Ri 



(A4) 



Namely, A R (r Q ) is the value of the random field at the position r a after applying a Gaussian kernel with a smoothing 
scale R a . The angular brackets denote an ensemble average used to calculate the autocorrelation matrix of the 
constraints and the cross correlation matrix of the constraints and the underlying 6 field. The constraints auto- 
correlation matrix is given by ^A Rq a r, \ = £a/3(l r a — r p\)y where the constraints auto-correlation function is 

= TSZui [d 3 kP(k)exp 



(2tt) s 



(A5) 



where P(k) is the power spectrum. Similarly, the constraints field cross-correlation matrix is given by 
£ a (|r — r |), where the constraints auto-correlation function is given by: 



<5(r) A 



1 



(2tt) s 



d 3 fc P{k) exp 



— ir ■ k — 



k 2 Ri 



Ra / = 



(A6) 



One should note that a constrained field is "effectively" constrained not only at the place where the constraints 
have been imposed, but also throughout the whole region within their respective correlation length. Regions located 
beyond that length obey the random part of the procedure. In principle, this can add other kind of structures that 
can affect the dynamical evolution of the imposed constraints. 

COMPARISON OF THE FTM CODE WITH THE SANTA BARBARA CLUSTER 



"The Santa Barbara cluster comparison project" l|Frenk et al.lll999T) was envisioned as a reliable comparison test 
between different cosmological numerical codes. It represents a hydrodynamical simulation of the formation of a rich 
cluster of galaxies. This exercise has become a standard test for new numerical codes in order to check for broad 
consistencies with other codes. Here we limit ourself to present only a comparison with respect to the DM properties 
of such cluster. The FTM has been tested extensively with respect to other dynamical codes. 

We have simulated the Santa Barbara cluster with a numerical r esol ution of 128 3 particles with a homogeneous 
sampling from the original 256 3 realization. The left-panel of Figure lA"Tl shows a density cut of t he final output tim e 
from the simulation. This figure has been constructed in the same way as those presented in iFrenk et al.l l)1999[) . 
i.e., , the image covers the inner 8 Mpc of the simulation cube, and it has been smoothed with a Gaussian kernel of 
i?G = 250 kpc. A visual comparison shows that our simulation possesses the same characteristics as the Santa Barbara 
one. This is further corroborated with the computed virial quantities: M v i r = 1-06 x 10 15 Mq, i? v ir = 2. 63 M pc, which 
are in good agreement with those reported in the Santa Barbara project. In the right panel of Fig. IA1I the radial 
density distribution of our cluster is presented. The red triangles represent our measurements and the red line shows 
the corresponding NFW profile. This distribution is in good agreement with respect to the averaged density profile 
reported in the Santa Barbara project. 
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